#include <iostream>
#include <conio.h>
#include <time.h>
#include <math.h>
#include <vector>

using namespace std;


/*
 * Generate a random number using rand() method
 * Return an odd number
 */
unsigned int GenerateRandomNumber(void)
{	
	return rand() | 1;
}

/*
 * Generate a random number using rand() method within n bits in decimal
 * Return an odd number
 */
unsigned int GenerateRandomNumber(int n)
{	
	int start = (int)pow(10.0, n -1);
	int end = start * 10;
	int u = (int)((double)rand() / (RAND_MAX + 1) * (end - start)
            + start);
	if (n == 1) {
		u = u == 1 ? u + 1 : u;
		return u;
	}
    return u | 1;
}

/*
 * Generate a random number using rand() method between start and end
 * Notice the end can't be reached, return an odd number
 */
unsigned int GenerateRandomNumber(int start, int end)
{	
	int u = (int)((double)rand() / (RAND_MAX + 1) * (end - start)
            + start);
    return u | 1;
}


/*
 * If n = a^b for integers a > 0 and b > 1, then return true, or return false.
 * The first step of the AKS algorithm
 */
bool VerifyPower(unsigned int n)
{
	unsigned int a = 2;
	for (; a <= floor(sqrt(n * 1.0)); ++a) {
		unsigned int temp = n;
		while (true) {			
			if (temp % a != 0)
				break;
			temp /= a;
			if (temp == 1)
				return true;
		}
	} 
	return false;	
}


/*
 * Get the Greatest Common Divisor
 * The iterative version
 */
unsigned int gcd(unsigned int a, unsigned int n)
{
	if (n == 0)
		return a;
	return gcd(n, a%n);
}

/*
 * Get the multiplicative order. Attention: the multiplicative order may not exists
 */
unsigned int GetMultiplicativeOrder(unsigned int a, unsigned int n)
{
	if (gcd(a, n) != 1)
		return -1;
	unsigned int k = 0, result = 1;
	while (true) {
		++k;
		result = (result * (a % n)) % n;
		if (result == 1)
			break;		
	}
	return k;
}

/*
 * Find the smallest r such that Or(n) > log(n) * log(n), with Or(n) being GetMultiplicativeOrder(n, r)
 */
unsigned int FindSmallestR(unsigned int n)
{
	unsigned int r = 2, k;
	while (true) {
		if (gcd(r, n) == 1) {
			k = GetMultiplicativeOrder(n, r);					
			if (k > (log10(n * 1.0) / log10(2.0)) * (log10(n * 1.0) / log10(2.0)))		
				return r;
		}		
		++r;
	}
}


/*
 * Euler's totient function
 */
unsigned int totient(unsigned int r)
{
	unsigned int count = 0;
	for (unsigned int i = 1; i <= r; ++i) {
		if (gcd(i, r) == 1)
			++count;
	}
	return count;
}


/*
 * Polynomial modulo
 * The future work is implementing it using Fourier (FFT algorithm)
 */
bool PolynomialMultiModulo(unsigned int n, unsigned int r, unsigned int end)
{
	int *polyArray = new int[r + 1];
	int *result = new int[r + 1];
	int *temp = new int[r + 1];
	bool equal;
	for (int a = 1; a <= end; a++) {
		memset(polyArray, 0, sizeof(int) * (r + 1));
		memset(result, 0, sizeof(int) * (r + 1));
		polyArray[0] = a;
		polyArray[1] = 1;

		result[0] = 1;

		int i, j, resultIndex = 1, polyIndex = 2;		
		int tempIndex = 0;
		int c = n;
		while (c != 0) {
			if (c % 2 == 1) {
				memset(temp, 0, sizeof(int) * (r + 1));			
				tempIndex = 0;
				for (i = 0; i < resultIndex; i++) {
					for (j = 0; j < polyIndex; j++) {
						int exp = (i + j) % r;
						int coef = (result[i] * polyArray[j]) % n;
						temp[exp] = (temp[exp] + coef) % n;

						if ((exp + 1) >= tempIndex) 
							tempIndex = exp + 1;
					}
				}

				int *inter = temp;
				temp = result;
				result = inter;

				resultIndex = tempIndex;
			}

			c /= 2;

			memset(temp, 0, sizeof(int) * (r + 1));
			tempIndex = 0;
			for (i = 0; i < polyIndex; i++) {
				for (j = 0; j < polyIndex; j++) {
					int exp = (i + j) % r;
					int coef = (polyArray[i] * polyArray[j]) % n;
					temp[exp] = (temp[exp] + coef) % n;
					if ((exp + 1) >= tempIndex) 
						tempIndex = exp + 1;
				}
			}
			int *inter = temp;
			temp = polyArray;
			polyArray = inter;

			polyIndex = tempIndex;
		}

		int exp1 = n % r;
		int coef1 = a % n;
		equal = true;

		if (result[exp1] == 1 && result[0] == coef1 && exp1 == (resultIndex - 1)) {		
			for (int idx = 1; idx < resultIndex; idx++) {				
				if (result[idx] != 0) {
					equal = false;
					break;
				}				
			}
		}

		if (!equal)
			break;
	} // end the outermost for loop
	delete []result;
	delete []polyArray;
	delete []temp;

	return equal;
}


int main(int argc, char **argv)
{
	int testNum = 1;
	srand((unsigned int)time(NULL));	
	int n = 0;
	bool isComposite = false;
	int bits;
	printf("Usage: Enter the bits of the prime in deciaml and 0 for exiting\n");
	cin >> bits;
	while (bits != 0) {
		isComposite = false;
		n = GenerateRandomNumber(bits);
		if (VerifyPower(n)) {			
			continue;
		} 
				
		unsigned int r = FindSmallestR(n);		
		
		for (unsigned int i = 2; i <= r; ++i) {
			if (gcd(i, n) > 1 && gcd(i, n) < n) {				
				isComposite = true;
				break;
			}
		}
		if (isComposite) 
			continue;

		if (n <= r) {
			printf("%d\n", n);	
			cin>>bits;
			continue;
		}
		
		unsigned int end = (int)floor(sqrt(totient(r) * 1.0) * log10(n * 1.0) / log10(2.0));

		if (!PolynomialMultiModulo(n, r, end))
			continue;
		
		printf("%d\n", n);
		cin>>bits;
	}
	return 0;
}